function Q=post_lhood(psi,theta,ests,inv_V)

% Unpack
L=length(psi)/4;
gamma=[psi(1:L),psi(2*L+1:3*L)];
kappa=[psi(L+1:2*L),psi(3*L+1:4*L)];

b_mu=exp(theta(1))/(1+exp(theta(1)));
b_alpha=theta(2);
b_beta=theta(3);
b_sigma=exp(theta(4));
b_gamma_mn=theta(5);
b_gamma_sd=exp(theta(6));
w_mu=exp(theta(7))/(1+exp(theta(7)));
w_alpha=theta(8);
w_beta=theta(9);
w_sigma=exp(theta(10));
w_gamma_mn=theta(11);
w_gamma_sd=exp(theta(12));
corr_gamma=(exp(theta(13))-1)/(exp(theta(13))+1);
corr_sigma=(exp(theta(14))-1)/(exp(theta(14))+1);

gamma_V=[b_gamma_sd^2,corr_gamma*b_gamma_sd*w_gamma_sd;corr_gamma*b_gamma_sd*w_gamma_sd,w_gamma_sd^2];
sigma_V=[b_sigma^2,corr_sigma*b_sigma*w_sigma;corr_sigma*b_sigma*w_sigma,w_sigma^2];
inv_gamma_V=inv(gamma_V); inv_sigma_V=inv(sigma_V);
gamma_m=[b_gamma_mn,w_gamma_mn];
kappa_m=[b_alpha,w_alpha]+repmat([b_beta,w_beta],L,1).*gamma;

% First part of likelihood (parameters)
Q=0;
for l=1:L
   Q=Q+(gamma(l,:)-gamma_m)*inv_gamma_V*(gamma(l,:)-gamma_m)'+(kappa(l,:)-kappa_m(l,:))*inv_sigma_V*(kappa(l,:)-kappa_m(l,:))'; %#ok
end

% Second part of likelihood (measurement)
sigma=exp(kappa);
pi=zeros(L,2);
b=zeros(L,2);
for l=1:L
    pi(l,1)=signalcdf(gamma(l,1),b_mu,sigma(l,1));
    pi(l,2)=signalcdf(gamma(l,2),w_mu,sigma(l,2));
    b(l,1)=1/(1+((1-b_mu)/b_mu)*exp((0.5-gamma(l,1))/sigma(l,1)^2));
    b(l,2)=1/(1+((1-w_mu)/w_mu)*exp((0.5-gamma(l,2))/sigma(l,2)^2));
end

type1=zeros(L,2);
type2=zeros(L,2);
delta=zeros(L,2);
for l=1:L
    type1(l,1)=normcdf(-0.5/sigma(l,1)-sigma(l,1)*log(((1-b(l,1))*b_mu)/(b(l,1)*(1-b_mu))));
    type1(l,2)=normcdf(-0.5/sigma(l,2)-sigma(l,2)*log(((1-b(l,2))*w_mu)/(b(l,2)*(1-w_mu))));
    type2(l,1)=normcdf(-0.5/sigma(l,1)+sigma(l,1)*log(((1-b(l,1))*b_mu)/(b(l,1)*(1-b_mu))));
    type2(l,2)=normcdf(-0.5/sigma(l,2)+sigma(l,2)*log(((1-b(l,2))*w_mu)/(b(l,2)*(1-w_mu))));
    delta(l,1)=type1(l,1)*b_mu/pi(l,1);
    delta(l,2)=type1(l,2)*w_mu/pi(l,2);
end

params=[pi;delta];
Q=Q+(ests(:)-params(:))'*inv_V*(ests(:)-params(:));



